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The  application  of  volume  rendering  techniques  to  the  display  of  surfaces  from  sampled  scalar 
functions  of  three  spatial  dimensions  is  explored.  Fitting  of  geometric  primitives  to  the  sampled 
data  is  not  required.  Images  arc  formed  by  directly  shading  each  sample  and  projecting  it  onto  the 
picture  plane.  Surface  shading  calculations  are  performed  at  every  voxel  with  local  gradient  vec¬ 
tors  serving  as  surface  normals.  In  a  separate  step,  surface  classification  operators  are  applied  to 
obtain  a  partial  opacity  for  every  voxel.  Operators  that  detect  isovalue  contour  surfaces  and  region 
boundary  surfaces  are  presented.  Independence  of  shading  and  classification  calculations  insures  an 
un distorted  visualization  of  3-D  shape.  Non-binary  classification  operators  insure  that  small  or 
poorly  defined  features  are  not  lost.  The  resulting  colors  and  opacities  are  composited  from  back  to 
front  along  viewing  rays  to  form  an  image.  The  technique  is  simple  and  fast,  yet  displays  surfaces 
exhibiting  smooth  silhouettes  and  few  other  aliasing  artifacts.  The  use  of  selective  blurring  and 
super-sampling  to  further  improve  image  quality  is  also  described.  Examples  from  two  applications 
are  given:  molecular  graphics  and  medical  imaging, 


1.  Introduction 

Visualization  of  scientific  computations  is  a  rapidly  growing  field  within  computer  graphics. 
A  large  subset  of  these  applications  involve  sampled  functions  of  three  spatial  dimensions,  also 
known  as  volume  data.  Surfaces  arc  commonly  used  to  visualize  volume  data  because  they  suc¬ 
cinctly  present  the  3-D  configuration  of  complicated  objects.  In  this  paper,  we  explore  the  use  of 
isovalue  contour  surfaces  to  visualize  electron  density  maps  for  molecular  graphics,  and  the  use  of 
region  boundary  surfaces  to  visualize  computed  tomography  (CT)  data  for  medical  imaging. 

The  currently  dominant  techniques  for  displaying  surfaces  from  volume  data  consist  of  apply¬ 
ing  a  surface  detector  to  the  sample  array,  fitting  geometric  primitives  to  the  detected  surfaces,  then 
rendering  these  primitives  using  conventional  surface  rendering  algorithms.  The  techniques  differ 
from  one  another  mainly  in  the  choice  of  primitives  and  the  scale  at  which  they  are  defined.  In  the 
medical  imaging  field,  a  common  approach  is  to  apply  thresholding  to  the  volume  data.  The  result¬ 
ing  binary  representation  can  be  rendered  by  treating  1-voxels  as  opaque  cubes  having  six  polygo¬ 
nal  faces  [1].  If  this  binary  representation  is  augmented  with  the  local  grayscale  gradient  at  each 
voxel,  substantial  improvements  in  surface  shading  can  be  obtained  [2-5].  Alternatively,  edge 
tracking  can  be  applied  on  each  slice  to  yield  a  set  of  contours  defining  features  erf  interest,  then  a 
mesh  of  polygons  can  be  constructed  connecting  the  contours  on  adjacent  slices  [6].  As  the  scale 
of  voxels  approaches  that  of  display  pixels,  it  becomes  feasible  to  apply  a  local  surface  detector  at 
each  sample  location.  This  yields  a  very  large  collection  of  voxel-sized  polygons,  which  can  be 
rendered  using  standard  algorithms  [7].  In  the  molecular  graphics  field,  methods  for  visualizing 
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electron  density  maps  include  stacks  of  isovalue  contour  lines,  ridge  lines  arranged  in  3-space  so  as 
to  connect  local  maxima  [8],  and  basket  meshes  representing  isovalue  contour  surfaces  [9]. 

All  of  these  techniques  suffer  from  the  common  problem  of  having  to  make  a  binary 
classification  decision:  either  a  surface  passes  through  the  current  voxel  or  it  does  not.  As  a  result, 
these  methods  often  exhibit  false  positives  (spurious  surfaces)  or  false  negatives  (erroneous  holes  in 
surfaces),  particularly  in  the  presence  of  small  or  poorly  defined  features. 

To  avoid  these  problems,  researchers  have  begun  exploring  the  notion  of  volume  rendering 
wherein  the  intermediate  geometric  representation  is  omiued.  Images  are  formed  by  shading  all 
data  samples  and  projecting  them  onto  the  picture  plane.  The  lack  of  explicit  geometry  does  not 
preclude  the  display  of  surfaces,  as  will  be  demonstrated  in  this  paper.  The  key  improvement 
offered  by  volume  rendering  is  that  it  provides  a  mechanism  for  displaying  weak  or  fuzzy  surfaces. 
This  capability  allows  us  to  relax  the  requirement,  inherent  when  using  geometric  representations, 
that  a  surface  be  either  present  or  absent  at  a  given  location.  This  in  turn  frees  us  from  the  neces¬ 
sity  of  making  binary  classification  decisions.  Another  advantage  of  volume  rendering  is  that  it 
allows  us  to  separate  shading  and  classification  operations.  This  separation  implies  that  the  accu¬ 
racy  of  surface  shading,  hence  the  apparent  orientation  of  surfaces,  does  not  depend  on  the  success 
or  failure  of  classification.  This  robustness  can  be  contrasted  with  rendering  techniques  in  which 
only  voxels  lying  on  detected  surfaces  are  shaded.  In  such  systems,  any  errors  in  classification 
result  in  incorrectly  oriented  surfaces. 

Smith  has  written  an  excellent  introduction  to  volume  rendering  [10],  Its  application  to  CT 
data  has  been  demonstrated  by  PIXAR  [11],  but  no  details  of  their  approach  have  been  published. 
The  technique  described  in  this  paper  grew  out  of  the  author’s  earlier  work  on  the  use  of  points  as 
a  rendering  primitive  [12].  Its  application  to  CT  data  was  first  reported  in  June,  1987  [13],  and  was 
presented  at  the  SP1E  Medical  Imaging  II  conference  in  January,  1988  [14]. 


2.  Rendering  pipeline 

The  volume  rendering  pipeline  used  in  this  paper  is  summarized  in  figure  1.  We  begin  with 
an  array  of  acquired  values  /0(x,)  at  voxel  locations  Xj  =  The  first  step  is  data  preparation 

which  may  include  correction  for  non-orthogonal  sampling  grids  in  electron  density  maps,  correc¬ 
tion  for  patient  motion  in  CT  data,  contrast  enhancement,  and  interpolation  of  additional  samples. 
The  output  of  this  step  is  an  array  of  prepared  values  /i(X|).  This  array  is  used  as  input  to  the  shad¬ 
ing  model  described  in  section  3,  yielding  an  array  of  voxel  colors  cx(xt),  X  =  rjjb.  In  a  separate 
step,  the  array  of  prepared  values  is  used  as  input  to  one  of  the  classification  procedures  described 
in  section  4,  yielding  an  array  of  voxel  opacities  a(xi).  Rays  are  then  cast  into  these  two  arrays 
from  the  observer  eyepoint.  For  each  ray,  a  vector  of  sample  colors  cx(X|)  and  opacities  o^q)  is 
computed  by  re- sampling  the  voxel  database  at  K  evenly  spaced  locations  X|  =  (x-o^Zf)  along  the 
ray  and  tri-linearly  interpolating  from  the  colors  and  opacities  in  the  eight  voxels  closest  to  each 
sample  location  as  shown  in  figure  2.  Finally,  a  fully  opaque  background  of  color  is  draped 
behind  die  dataset  and  the  re-sampled  colors  and  opacities  are  merged  with  each  other  and  with  the 
background  by  compositing  in  back-to-front  order  to  yield  a  single  color  Cx(uj)  for  the  ray,  and, 
since  only  one  ray  is  cast  per  image  pixel,  for  the  pixel  location  14  a  (u:,vj)  as  well. 

The  compositing  calculations  referred  to  above  are  simply  linear  interpolations.  Specifically, 
the  color  CMlx(ti|)  of  the  ray  as  it  leaves  each  sample  location  is  related  to  the  color  C^uj)  of  the 
ray  as  it  enters  and  the  color  cx(xj)  and  opacity  a (*j)  at  that  sample  location  by  the  transparency 
formula 


l)  *  C*a(“|)(l  -  afx]))  +  c^x^aty). 

Solving  for  pixel  color  Cx(uj)  in  terms  of  the  vector  of  sample  colors  cx(*j)  and  opacities  a(x|) 
along  the  associated  viewing  ray  gives 
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C\(uj)  =  Cx(u;,vp  = 
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where  cx(x;0';^o)  =  e**,.*  and  a(x;,yj,z0)  =  1. 


(1) 


3.  Shading 

Using  the  rendering  pipeline  presented  above,  the  mapping  from  acquired  data  to  color  pro¬ 
vides  3-D  shape  cues,  but  does  not  participate  in  the  classification  operation.  Accordingly,  a  shad¬ 
ing  model  was  selected  that  provides  a  satisfactory  illusion  of  smooth  surfaces  at  a  reasonable  cost 
It  is  not  the  main  point  of  the  paper  and  is  presented  mainly  for  completeness.  The  model  chosen 
is  due  to  Phong  [IS]: 

Cx(x  l)  =  + 

^^10‘WWHr]  (2) 

where 

cx(X|)  =  X’th  component  of  color  at  voxel  location  x,,  X  =  r,gjb, 
cPi  =  X’th  component  of  color  of  parallel  light  source, 

=  ambient  reflection  coefficient  for  X’th  color  component, 
k<a  =  diffuse  reflection  coefficient  for  X’th  color  component, 
kt  x  =  specular  reflection  coefficient  for  X’th  color  component. 


n  =  exponent  used  to  approximate  highlight, 
kx,  k2  =  constants  used  in  linear  approximation  of  depth-cueing, 
dfxO  =  perpendicular  distance  from  picture  plane  to  voxel  location  X|, 
N(xi)  =  surface  normal  at  voxel  location  X|, 

L  =  normalized  vector  in  direction  of  light  source, 

H  =  normalized  vector  in  direction  of  maximum  highlight. 

Since  a  parallel  light  source  is  used,  L  is  a  constant.  Furthermore, 


H  = 


V  +  L 
IV  +  LI 


where 
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V  =  normalized  vector  in  direction  of  observer. 

Since  an  orthographic  projection  is  used,  V  and  hence  H  are  also  constants.  Finally,  the  surface 
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normal  is  given  by 


N(xi)  = 


V/frt) 


where  the  gradient  vector  V/(xq  is  approximated  using  the  operator 

VJl*d  =  = 


4.  Classification 

The  mapping  from  acquired  data  to  opacity  performs  the  essential  task  of  surface 
classification.  We  will  first  consider  the  rendering  of  isovalue  contour  surfaces  in  electron  density 
maps,  i.e.  surfaces  defined  by  points  of  equal  electron  density.  We  will  then  consider  the  render- 
ing  of  region  boundary  surfaces  in  computed  tomography  (CT)  data,  i.e.  surfaces  bounding  tissues 
of  constant  CT  number. 


4.1.  Isovalue  contour  surfaces 

Determining  the  structure  of  large  molecules  is  a  difficult  problem.  The  method  most  com¬ 
monly  used  is  ab  initio  interpretation  of  electron  density  maps,  which  represent  the  averaged  den¬ 
sity  of  a  molecule's  electrons  as  a  function  of  position  in  3-space.  These  maps  are  obtained  from 
X-ray  diffraction  studies  of  crystallized  samples  of  the  molecule. 

One  obvious  way  to  display  isovalue  surfaces  is  to  opaquely  render  all  voxels  having  values 
greater  than  some  threshold.  This  produces  3-D  regions  of  opaque  voxels  the  outermost  layer  of 
which  is  the  desired  isovalue  surface.  Unfortunately,  this  solution  prevents  display  of  multiple  con¬ 
centric  semi-transparent  surfaces,  a  very  useful  capability.  Using  a  window  in  place  of  a  threshold 
does  not  solve  the  problem.  If  die  window  is  too  narrow,  holes  appear.  If  it  too  wide,  display  of 
multiple  surfaces  is  constrained.  In  addition,  the  use  of  thresholds  and  windows  introduces  artifacts 
into  the  image  that  are  not  present  in  the  data. 

The  classification  procedure  employed  in  this  study  begins  by  assigning  an  opacity  a,  to  vox¬ 
els  having  selected  value  /„  and  assigning  an  opacity  of  zero  to  aU  other  voxels.  In  order  to  avoid 
aliasing  artifacts,  we  would  also  like  voxels  having  values  close  to  /,  to  be  assigned  opacities  close 
to  a,.  The  most  pleasing  image  is  obtained  if  the  thickness  of  this  transition  region  stays  constant 
throughout  the  volume.  We  approximate  this  effect  by  having  the  opacity  fall  off  as  we  move 
away  from  the  selected  value  at  a  rate  inversely  proportional  to  the  magnitude  of  the  local  gradient 
vector. 

This  mapping  is  implemented  using  the  expression 
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a(X|)  =  a,  < 


r  |  I 
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if  IV/(x,)l  =  0  and 

/W  =/. 
if  IV/(xi)l  >  0  and 
f[xd  -  r  IV/x,)!  S/,  £ 
/W  +  r  IVflx,)! 
otherwise 


(3) 


where  r  is  the  desired  thickness  in  voxels  of  the  transition  region,  and  the  gradient  vector  is  approx¬ 
imated  using  the  operator  given  in  section  3.  A  graph  of  a(xi)  as  a  function  of  /(x,)  and  for 
typical  values  of /„  a*,  and  r  is  shown  in  figure  3. 

If  more  than  one  isovalue  surface  is  to  be  displayed  in  a  single  image,  they  can  be  classified 
separately  and  their  opacities  combined.  Specifically,  given  selected  values 
/v  n  =  1, ...  .N,  N  £  1,  opacities  a,a  and  transition  region  thicknesses  r„  we  can  use  equation  (3) 
to  compute  an(xl),  then  apply  the  relation 

a pM)  =  1  -  FI  (1  -  a,(Xi)). 


42.  Region  boundary  surfaces 

From  a  densitometric  point  of  view,  the  human  body  is  a  complex  arrangement  of  biological 
tissues  each  of  which  is  fairly  homogeneous  and  of  predictable  density.  Clinicians  are  mostly 
interested  in  the  boundaries  between  tissues,  from  which  the  sizes  and  spatial  relationships  of  ana¬ 
tomical  features  can  be  inferred. 

Although  many  researchers  use  isovalue  surfaces  for  the  display  of  medical  data,  it  is  not 
clear  that  they  are  well  suited  for  that  purpose.  The  reason  can  be  explained  briefly  as  follows. 
Given  an  anatomical  scene  containing  two  tissue  types  A  and  B  having  values  f,A  and  where 
At  < /**•  d*1*  acquisition  will  produce  voxels  having  values  such  that  /„  S  fy$.  Thin 

features  of  tissue  type  B  may  be  represented  by  regions  in  which  all  voxels  bear  values  less  than 
/„ .  Indeed,  there  is  no  threshold  value  greater  than  guaranteed  to  detect  arbitrarily  thin  regions 

of  type  B,  and  thresholds  close  to  /„  are  as  likely  to  detect  noise  as  signal. 

The  procedure  employed  in  this  study  is  based  on  the  following  simplified  model  of  anatomi¬ 
cal  scenes  and  the  CT  scanning  process.  We  assume  that  scenes  contain  an  arbitrary  number  of  tis¬ 
sue  types  bearing  CT  numbers  falling  within  a  small  neighborhood  of  some  known  value.  We 
further  assume  that  tissues  of  each  type  touch  tissues  of  at  most  two  other  types  in  a  given  scene. 
Finally,  we  assume  that,  if  we  order  the  types  by  CT  number,  then  each  type  touches  only  types 
adjacent  to  it  in  the  ordering.  Formally,  given  N  tissue  types  bearing  CT  numbers 
/, ,  n  *  1, . . .  Jf,  N  2  1  such  that  /,  </„  ,,  m  »  1, . . .  Ji-\,  then  no  tissue  of  CT  number  /„ 

a  a  »j 

touches  any  tissue  of  CT  number/,^,  >  1. 

If  these  criteria  are  met,  each  tissue  type  can  be  assigned  an  opacity  and  a  piecewise  linear 
mapping  can  be  constructed  that  convent  voxel  value  /%  to  opacity  aY<1  voxel  value /V|  to  opacity 
a,^,  and  intermediate  voxel  values  to  intermediate  opacities.  Note  that  all  voxels  are  typically 
mapped  to  some  non- zero  opacity  and  will  thus  contribute  to  the  final  image.  This  scheme  insures 
that  thin  regions  of  tissue  will  still  appear  in  the  image,  even  if  only  as  faint  wisps.  Note  also  that 
violation  of  die  adjacency  criteria  leads  to  voxels  that  cannot  be  unambiguously  classified  as 
belonging  to  one  region  boundary  or  another  and  hence  cannot  be  rendered  correctly  using  this 
method. 
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The  superimposition  of  multiple  semi-transparent  surfaces  such  as  skin  and  bone  can  substan¬ 
tially  enhance  the  comprehension  of  CT  data.  In  order  to  obtain  such  effects  using  volume  render¬ 
ing,  we  would  like  to  suppress  the  opacity  of  tissue  interiors  while  enhancing  the  opacity  of  their 
bounding  surfaces.  We  implement  this  by  scaling  the  opacities  computed  above  by  the  magnitude 
of  the  local  gradient  vector. 

Combining  these  two  operations,  we  obtain  a  set  of  expressions 


a(x.)  =  IV/WI 


m-L 


/%..  "A 


if  A.  Sfa) 


otherwise 


(4) 


for  n  *  1, .  .  .  /M,  N  £  1.  The  gradient  vector  is  approximated  using  the  operator  given  in  sec¬ 
tion  3.  A  graph  of  a(x()  as  a  function  of  f[x{)  and  iV^x^l  for  three  tissue  types  A,  B,  and  C,  having 
typical  values  of  /v/„#,/,c,  OvA.  a^,  and  otvc  is  shown  in  figure  4. 


5.  Discussion 


5.1.  Computational  complexity 

One  of  the  strengths  of  the  rendering  method  presented  in  this  paper  is  its  modularity.  By 
storing  intermediate  results  at  various  places  in  the  pipeline,  the  cost  of  generating  a  new  image 
after  a  change  in  shading  or  classification  parameters  is  minimized.  Let  us  consider  some  typical 
cases. 

Given  acquired  value  /(xi)  and  gradient  magnitude  IVftxi)l,  the  mapping  to  opacity  a(xi)  can 
be  implemented  with  one  lookup  table  reference.  This  implies  that  if  we  store  gradient  magnitudes 
for  all  voxels,  computation  of  new  opacities  following  a  change  in  classification  parameters  entails 
only  generation  of  a  new  lookup  table  followed  by  one  table  reference  per  voxel. 

The  cost  of  computing  new  colors  c*(x, )  following  a  change  in  observer  direction  V,  light 
source  direction  L,  or  other  shading  parameter  is  more  substantial.  Effective  rotation  sequences 
can  be  generated,  however,  using  a  single  set  of  colors.  The  visual  manifestation  of  fixing  the 
shading  is  that  light  sources  appear  to  travel  around  with  the  data  as  it  routes  and  highlights  are 
incorrect  Since  we  are  visualizing  imaginary  or  invisible  phenomena  anyway,  observers  are  sel¬ 
dom  troubled  by  this  effect 

The  most  efficient  way  to  produce  a  rotation  sequence  is  then  to  hold  both  colors  and  opaci¬ 
ties  constant  and  alter  only  the  direction  in  which  rays  are  cast  If  we  assume  a  square  image  n 
pixels  wide  and  use  orthographic  projection,  in  which  case  sample  coordinates  can  be  efficiently 
calculated  using  differencing,  the  combined  cost  of  ray  tracing,  re-sampling,  and  compositing  to 
compute  n 2  pixels  is  IKn2  additions,  2 Kn2  tri-linear  interpolations,  and  Kir  linear  interpolations, 
where  K  is  the  number  of  sample  locations  along  each  ray. 


7 


5.2.  Image  quality 

Although  the  notation  used  in  equation  (1)  has  been  borrowed  from  the  literature  of  image 
compositing  [16],  the  analogy  is  not  exact,  and  the  differences  are  fundamental.  Volume  data  con¬ 
sists  of  samples  taken  from  a  bandlimited  3-D  scene,  whereas  the  data  acquired  from  an  image 
digitizer  consists  of  samples  taken  from  a  bandlimited  2-D  projection  of  a  3-D  scene.  Unless  we 
reconstruct  the  3-D  scene  that  gave  rise  to  our  volume  data,  we  cannot  compute  an  accurate  projec¬ 
tion  of  iL  Volume  rendering  performs  no  such  reconstruction.  Image  quality  is  therefore  limited 
by  the  number  of  viewing  rays.  In  the  current  implementation,  we  cast  one  ray  per  pixel.  Such 
point  sampling  would  normally  produce  strong  aliasing,  but,  by  using  non-binary  classification  deci¬ 
sions,  we  carry  much  of  the  bandlimiting  inherent  in  the  acquired  data  over  into  the  image,  sub¬ 
stantially  reducing  aliasing  artifacts.  Stated  another  way,  we  are  depending  on  3-D  bandlimiting  to 
avoid  aliasing  in  2-D  projections. 

Within  these  limitations,  there  are  two  ways  to  improve  image  quality,  blurring  and  super- 
sampling.  If  the  array  of  acquired  values  are  bluned  slightly  during  data  preparation,  the  oversharp 
surface  silhouettes  occasionally  exhibited  by  volume  renderings  are  softened.  Alternatively,  we  can 
apply  blurring  to  the  opacities  generated  by  the  classification  procedure,  but  leave  the  shading 
untouched.  This  has  the  effect  of  softening  silhouettes  without  adversely  affecting  the  crispness  of 
surface  detail. 

The  decision  to  reduce  aliasing  at  the  expense  of  resolution  arises  from  two  conflicting  goals: 
generating  artifact-free  images  and  keeping  rendering  costs  low.  In  practice,  the  slight  loss  in 
image  sharpness  might  not  be  disadvantageous.  Indeed,  it  is  not  clear  that  the  accuracy  afforded  by 
more  expensive  visibility  calculations  is  useful,  at  least  for  the  types  of  data  considered  in  this 
study.  Blurry  silhouettes  have  less  visual  impact,  but  they  reflect  the  true  imprecision  in  our 
knowledge  of  surface  locations. 

An  alternative  means  for  improving  im^ge  quality  is  super-sampling.  The  basic  idea  is  to 
interpolate  additional  samples  between  the  acquired  ones  prior  to  compositing.  If  the  interpolation 
method  is  a  good  one,  the  accuracy  of  the  visibility  calculations  is  improved,  reducing  some  kinds 
of  aliasing.  Another  option  is  to  apply  this  interpolation  during  data  preparation.  Although  this 
alternative  substantially  increases  computational  expense  throughout  the  remainder  of  the  pipeline, 
it  improves  the  accuracy  of  our  shading  and  classification  calculations  as  well  as  our  visibility. 


6.  Implementation  and  results 

The  dataset  used  in  the  molecular  graphics  study  is  a  113  x  113  x  113  voxel  portion  of  an 
electron  density  map  for  the  protein  Cytochrome  B5.  Figure  S  shows  four  slices  spaced  10  voxels 
apart  in  this  dataset.  Each  whitish  cloud  represents  a  single  atom.  Using  the  shading  and 
classification  calculations  described  in  sections  3  and  4.1,  colors  and  opacities  were  computed  for 
each  voxel  in  the  expanded  dataset  These  calculations  required  5  minutes  on  a  SUN  4/280  having 
32MB  of  main  memory.  Ray  tracing  and  compositing  were  performed  as  described  in  section  2 
and  took  40  minutes,  yielding  the  image  in  figure  6. 

The  dataset  used  in  the  medical  imaging  study  is  a  CT  study  of  a  cadaver  and  was  acquired 
as  113  slices  of  256  x  256  samples  each.  Using  the  shading  and  classification  calculations 
described  in  sections  3  and  4.2,  two  sets  of  colon  and  opacities  were  computed,  one  showing  the 
air-skin  interface  and  a  second  showing  the  tissue-bone  interface.  The  computation  of  each  set 
required  20  minutes.  Using  the  compositing  algorithm  described  in  section  2,  two  views  were  then 
computed  from  each  set  of  colon  and  opacities,  producing  four  images  in  all,  as  shown  in  figure  7. 
The  compulation  of  each  view  required  an  additional  20  minutes.  The  horizontal  bands  through  the 
patient's  teeth  in  aU  of  these  images  are  aitifrcts  due  to  scattering  of  X-rays  from  dental  fillings 
and  are  preaent  in  the  acquired  data.  The  bonds  across  her  forehead  and  under  her  chin  in  the  air- 


skin  images  are  gauze  bandages  used  to  immobilize  her  head  during  scanning.  Her  skin  and  nose 
cartilage  are  rendered  semi-transparently  over  the  bone  surface  in  the  tissue-bone  images. 


Figure  8  was  generated  by  combining  halves  from  each  of  the  two  sets  of  colots  and  opaci¬ 
ties  already  computed  for  figure  7.  Heightened  transparency  of  the  temporal  bone  and  the  bones 
surrounding  the  maxillary  sinuses  -  more  evident  in  moving  sequences  than  in  a  static  view  •  is  due 
to  generalized  osteoporosis.  It  is  worth  noting  that  rendering  techniques  employing  binary 
classification  decisions  would  likely  display  holes  here  instead  of  thin,  wispy  surfaces. 

The  dataset  used  in  figures  9  and  10  is  of  the  same  cadaver,  but  was  acquired  as  1 13  slices  of 
512  x  512  samples  each.  Figure  9  was  generated  using  the  same  procedure  as  for  figure  7,  but 
casting  four  rays  per  slice  in  the  vertical  direction  in  order  to  correct  for  the  aspect  ratio  of  ihe 
dataset.  Figure  10  was  generated  by  expanding  the  dataset  to  452  slices  using  a  cubic  B-spline  in 
the  vertical  direction,  then  generating  an  image  from  the  larger  dataset  by  casting  one  ray  per  slice. 
As  expected,  more  detail  is  apparent  in  figure  10  than  figure  9. 


7.  Conclusions 

Volume  rendering  has  been  shown  to  be  an  effective  modality  for  the  display  of  surfaces  from  sam¬ 
pled  scalar  functions  of  three  spatial  dimensions.  As  demonstrated  by  the  figures,  it  can  generate 
images  exhibiting  approximately  equivalent  resolution,  yet  fewer  interpretation  errors,  than  tech¬ 
niques  relying  on  geometric  primitives. 

Despite  its  advantages,  volume  rendering  has  several  problems.  The  omission  of  an  inter¬ 
mediate  geometric  representation  makes  selection  of  appropriate  shading  parameters  critical  to  the 
effectiveness  of  the  visualization.  Slight  changes  in  opacity  ramps  or  interpolation  methods  radi¬ 
cally  alter  the  features  that  are  seen  as  well  as  the  overall  quality  of  the  image.  For  example,  the 
thickness  of  the  transition  region  surrounding  the  isovalue  contour  surfaces  described  in  section  4.1 
stays  constant  only  if  the  local  gradient  magnitude  stays  constant  within  a  radius  of  r  voxels  around 
each  point  on  the  surface.  The  time  and  ensemble  averaging  inherent  in  X-ray  crystallography  usu¬ 
ally  yields  suitable  data,  but  there  axe  considerable  variations  among  datasets.  Algorithms  are 
needed  that  automatically  select  an  optimum  value  for  r  based  on  the  characteristics  of  a  particular 
dataset. 

Volume  rendering  is  also  very  sensitive  to  artifacts  in  the  acquisition  process.  For  example, 
CT  scanners  generally  have  anisotropic  spatial  sensitivity.  This  problem  manifests  itself  as  striping 
in  images.  With  live  subjects,  patient  motion  is  also  a  serious  problem.  Since  shading  calculations 
are  strongly  dependent  on  the  orientation  of  the  local  gradient,  slight  mis-alignmems  between  adja¬ 
cent  slices  produce  strong  striping. 

An  issue  related  specifically  to  region  boundary  surfaces  is  the  rendering  of  datasets  not 
meeting  the  adjacency  criteria  described  in  section  4.2.  This  includes  most  internal  soft  tissue 
organs  in  CT  data.  One  simple  solution  would  be  for  users  to  interactively  clip  or  carve  away 
unwanted  tissues,  thus  isolating  subsets  of  the  acquired  data  that  meet  the  adjacency  criteria.  Since 
the  user  or  algorithm  is  not  called  upon  to  define  geometry,  but  merely  to  isolate  regions  of 
interest,  this  approach  promises  to  be  easier  and  to  produce  better  images  than  techniques  involving 
surface  fitting.  A  more  sophisticated  approach  would  be  to  combine  volume  rendering  with  high- 
level  object  definition  methods  in  an  interactive  setting.  Initial  visualizations,  made  without  the 
benefit  of  object  definition,  would  be  used  to  guide  scene  analysis  snd  segmentation  algorithms, 
which  would  in  turn  be  used  to  isolate  regions  of  interest,  producing  s  better  visualization.  If  the 
output  of  such  segmentation  algorithms  included  confidence  levels  or  probabilities,  they  could  be 
mapped  to  opacity  and  thus  modulate  the  appearance  of  the  image. 

Although  this  paper  focuses  on  display  of  surfaces,  the  same  pipeline  cm  be  easily  modifed 
to  render  interstitial  volumes  as  semi-transparent  gels.  Color  snd  texture  cm  be  added  to  represent 
such  variables  as  gradient  magnitude  or  araaomkal  tissue  type.  Visualizations  combining  sampled 
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and  geometric  data  also  hold  much  promise.  For  example,  it  might  be  useful  to  superimpose  ball- 
and-stick  molecular  models  onto  electron  density  maps  or  medical  prostheses  onto  CT  scans.  To 
obtain  correct  visibility,  a  true  3-D  merge  of  the  sampled  and  geometric  data  must  be  performed. 
One  possible  solution  is  to  scan-convert  the  geometry  directly  into  the  sampled  array  and  render  the 
ensemble.  Alternatively,  classical  ray  tracing  of  the  geometry  can  be  incorporated  direcdy  into  the 
volume  rendering  pipeline.  Another  useful  tool  would  be  the  ability  to  perform  a  true  3-D  merge 
of  two  or  more  visualizations,  allowing,  for  example,  the  superimposition  of  radiation  treatment 
planning  isodose  surfaces  over  CT  data. 

The  prospects  for  rapid  interactive  manipulation  of  volume  data  are  encouraging.  By  pre- 
computing  colors  and  opacities  and  storing  them  in  intermediate  3-D  datasets,  we  simplify  the 
image  generation  problem  to  one  of  geometrically  transforming  two  values  per  voxel  and  composit¬ 
ing  the  results.  One  promising  technique  for  speeding  up  these  transformations  is  to  combine  a  3- 
pass  version  of  2-pass  texture  mapping  [17]  with  compositing.  By  re-sampling  separately  in  each 
of  three  orthogonal  directions,  computational  expense  is  reduced.  Given  a  suitably  interpolated 
sample  array,  it  might  be  possible  to  omit  the  third  pass  entirely,  compositing  transformed  2-D 
slices  together  to  form  an  image.  This  further  suggests  that  hardware  implementations  might  be 
feasible.  A  recent  survey  of  architectures  for  rendering  voxel  data  is  given  by  Kaufman  [18].  A 
suitably  interconnected  2-D  array  of  processors  with  sufficient  backing  storage  might  be  capable  of 
producing  visualizations  of  volume  datasets  in  real-time  or  near  real-time. 
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Figure  1  -  Overview  of  volume  rendering  pipeline 


Figure  2  -  Ray  tracing  /  re-sampling  steps 


opacity  a(Xj) 


Figure  3  -  Calculation  of  opacities  for  isovalue  contour  surfaces 


opacity  a(x.) 


Figure  4  -  Calculation  of  opacities  for  region  boundary  surfaces 
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